{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<img src=\"../common/rfsoc_book_banner.jpg\" alt=\"University of Strathclyde\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Notebook Set H\n",
    "\n",
    "---\n",
    "\n",
    "## 03 - FEC Channel Simulation\n",
    "This is the third notebook in the series exploring Soft Decision Forward Error Correction (SD-FEC) on RFSoC. In previous notebooks we have learned the underlying principals of forward error correction and have used an SD-FEC integrated block to encode data in the programmable logic of the RFSoC. In this notebook we will simulate an additive white Gaussian noise (AWGN) channel for our encoded data to pass through. \n",
    "\n",
    "## Table of Contents\n",
    "* [1. Introduction](#nb3_introduction)\n",
    "    * [1.1. Overview](#nb3_overview)\n",
    "    * [1.2. Notebook Setup](#nb3_notebook_setup)\n",
    "* [2. Symbol Mapping (Baseband Modulation)](#symbol_mapping)\n",
    "* [3. AWGN](#awgn)\n",
    "    * [3.1. Calculate Noise Variance](#calculate_noise_variance)\n",
    "    * [3.2. Generate Noise Signal](#generate_noise_signal)\n",
    "    * [3.3. Add Noise](#add_noise)\n",
    "* [4. Log Likelihood Ratio (Soft Demodulation)](#llr)\n",
    "* [5. Conclusion](#nb3_conclusion)\n",
    "\n",
    "## References\n",
    "* [1] - [AMD-Xilinx, \"Soft-Decision FEC Integrated Block v1.1: LogiCORE IP Product Guide\", October 2022](https://docs.xilinx.com/r/en-US/pg256-sdfec-integrated-block)\n",
    "* [2] - [Krishna Sankar, \"Softbit for 16QAM\", July 2009](http://www.dsplog.com/2009/07/05/softbit-16qam/)\n",
    "\n",
    "## Revision\n",
    "* **v1.0** | 17/01/23 | *First Revision*\n",
    "\n",
    "---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction <a class=\"anchor\" id=\"nb3_introduction\"></a>\n",
    "The purpose of error correction codes is to identify and fix errors in a data block that result from damage or noise. As these notebooks are focused on radio applications using RFSoC, our encoded data will be exposed to an AWGN channel as part of a simplified radio pipeline. Our encoded data will first be baseband modulated using a 16-QAM scheme, then exposed to an AWGN channel, and finally baseband demodulated using soft decisions. The resulting data will be noisy encoded data in the form of soft bits.\n",
    "\n",
    "### 1.1. Overview <a class=\"anchor\" id=\"nb3_overview\"></a>\n",
    "Figure 1 illustrates how this notebook is intended to interact with notebooks two and four. Encoded data from notebook two is subjected to our simplified radio pipeline in this notebook and the resulting soft bits or Log Likelihood Ratios (LLRs) are used to perform decoding in notebook four. Note that, while this notebook is intended to be executed as part of a larger group of notebooks on RFSoC, if you are running this on a personal computer, then an encoded data block will instead be imported from a file.\n",
    "\n",
    "<a class=\"anchor\" id=\"fig-1\"></a>\n",
    "<center><figure>\n",
    "<img src='./images/nb_3_overview.svg' width='500'/>\n",
    "    <figcaption><b>Figure 1: Functional block diagram illustrating the channel simulation pipeline.</b></figcaption>\n",
    "</figure></center>\n",
    "\n",
    "\n",
    "### 1.2. Notebook Setup <a class=\"anchor\" id=\"nb3_notebook_setup\"></a>\n",
    "We setup the notebook by importing the required Python libraries and the encoded data from the previous notebook (or file)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import strath_sdfec.helper_functions as hf\n",
    "\n",
    "%store -r tx_enc_buf rx_bits\n",
    "\n",
    "if 'rx_bits' not in locals():\n",
    "    print('No stored data found, importing from file..')\n",
    "    tx_enc_buf = np.loadtxt('data/docsis_init_ranging_tx.txt')\n",
    "    rx_bits = np.loadtxt('data/docsis_init_ranging_rx.txt')\n",
    "    \n",
    "hf.plot_samples('Encoded Serial Buffer',\n",
    "                [range(len(rx_bits))],\n",
    "                [rx_bits])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Baseband Modulation <a class=\"anchor\" id=\"symbol_mapping\"></a>\n",
    "The first stage in the simplified radio pipeline is baseband modulation. For digital data to be sent over a wireless communications channel, the serial bits must first be converted to voltage amplitude levels. This conversion is called baseband modulation or bit-to-symbol mapping. This process involves dividing the serial bits into groups of a certain length. These groups are referred to as symbols and each symbol will have an associated voltage level or levels in the case of complex signals. How the symbols are mapped to a voltage level is determined by the modulation scheme employed. In this notebook, we use the 16-QAM scheme and Gray coding, where adjacent symbols differ by only 1 bit. Figure 2 depicts a constellation of the symbol mapping used.\n",
    "\n",
    "<a class=\"anchor\" id=\"fig-2\"></a>\n",
    "<center><figure>\n",
    "<img src='./images/symbol_mapping.png' width='800'/>\n",
    "    <figcaption><b>Figure 2: Gray coded symbol mapping constellation for 16-QAM.</b></figcaption>\n",
    "</figure></center>\n",
    "\n",
    "By grouping bits together and assigning them a specific amplitude level, more data is able to be transmitted per sample. The more bits that form a symbol, the higher the transmission data rate is. However, you will notice that this also means that there is less space between symbols in the constellation, meaning the signal is more susceptible to noise. The 16-QAM scheme has a total of 16 symbols meaning each symbol is made up of 4 bits. From this constellation, we can easily see that the two most significant bits are responsible for the real amplitude level and the two least significant bits determine the imaginary amplitude level. \n",
    "\n",
    "We can write a function that loops through a serial bit-stream in groups of 4 bits and ascertains the real and imaginary amplitude levels by comparing the two MSBs and two LSBs to certain bit patterns as detailed in the constellation diagram. In the cell below, we do exactly that, and the function returns a baseband modulated complex signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def symbol_map(serial_data):\n",
    "    mod_I = np.empty(0)\n",
    "    mod_Q = np.empty(0)\n",
    "    \n",
    "    # Iterate through serial bits in increments of 4\n",
    "    for i in range(0, len(serial_data), 4):\n",
    "        symbol = serial_data[i : i+4]    # Group 4 bits\n",
    "\n",
    "        b3b2 = symbol[0:2]    # MSBs\n",
    "        b1b0 = symbol[2:4]    # LSBs\n",
    "\n",
    "        # Real amplitude\n",
    "        if np.array_equal(b3b2,[0,0]):\n",
    "            I = -3\n",
    "        elif np.array_equal(b3b2,[0,1]):\n",
    "            I = -1\n",
    "        elif np.array_equal(b3b2,[1,1]):\n",
    "            I = 1\n",
    "        elif np.array_equal(b3b2,[1,0]):\n",
    "            I = 3\n",
    "\n",
    "        # Imaginary amplitude\n",
    "        if np.array_equal(b1b0,[0,0]):\n",
    "            Q = 3\n",
    "        elif np.array_equal(b1b0,[0,1]):\n",
    "            Q = 1\n",
    "        elif np.array_equal(b1b0,[1,1]):\n",
    "            Q = -1\n",
    "        elif np.array_equal(b1b0,[1,0]):\n",
    "            Q = -3\n",
    "\n",
    "        mod_I = np.append(mod_I, I)\n",
    "        mod_Q = np.append(mod_Q, Q)\n",
    "        \n",
    "    signal = mod_I + 1j*mod_Q   \n",
    "    return signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the function defined above, we can plot the samples of our new signal as well as a constellation diagram. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = symbol_map(rx_bits)\n",
    "\n",
    "hf.plot_samples('Encoded Signal',\n",
    "                [range(len(signal)),range(len(signal))],\n",
    "                [np.real(signal),np.imag(signal)])\n",
    "\n",
    "hf.plot_constellation(signal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our encoded data has been converted from a stream of individual bits (0s and 1s) of length *L* to a complex 16-QAM signal of length *L*/4 that has 4 amplitude levels in both the real and imaginary domain. \n",
    "\n",
    "## 3. AWGN Channel <a class=\"anchor\" id=\"awgn\"></a>\n",
    "We will now subject our signal to noise that we can control the level of. To achieve this, we will create a second signal of normally distributed random numbers that will act as our noise and be added to our data signal. We must be able to adjust the spread of this new signal's amplitude such that, when added to our data signal, will result in a desired Signal-to-Noise Ratio (SNR). The following subsections outline this process.\n",
    "\n",
    "### 3.1. Calculate Noise Variance <a class=\"anchor\" id=\"calculate_noise_variance\"></a>\n",
    "Variance is the measure of spread around the mean. If we take the equation of Signal-to-Noise Ratio expressed in decibels:\n",
    "\n",
    "$$\n",
    "    \\mathit{SNR\\;(dB)} =  10\\log_{10} \\left(\\frac{\\sigma^2_{signal}}{\\sigma^2_{noise}}\\right)\n",
    "$$\n",
    "We can rearrange for noise variance:\n",
    "$$\n",
    "    \\sigma^2_{noise} = \\frac{\\sigma^2_{signal}}{10^{\\frac{\\mathit{\\;SNR\\;(dB)}}{10}}}\n",
    "$$\n",
    "\n",
    "This allows us to update the SNR value and obtain the variance that the noise must have to result in our desired SNR. The cell below performs this calculation. While the signal variance is measured here, it is known that for the 16-QAM scheme employed and the amplitude levels used, the resulting signal variance will be 10 provided the mapped bits are randomly generated with uniform distribution. The more bits measured over, the more accurate this result will be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNR = 15\n",
    "\n",
    "snr_ratio = 10**(SNR/10)\n",
    "var_signal = np.var(signal)\n",
    "var_noise = var_signal / snr_ratio\n",
    "\n",
    "print('Desired SNR: %d dB' % SNR)\n",
    "print('Signal Variance:', var_signal)\n",
    "print('Noise Variance:' ,var_noise)\n",
    "print('SNR (ratio):', var_signal/var_noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. Generate Noise Signal <a class=\"anchor\" id=\"generate_noise_signal\"></a>\n",
    "Our noise signal is to be complex and therefore the calculated variance applies to this complex signal. To create a complex noise signal with the calculated variance we must create two separate noise signals, one to act as the real component and the other the imaginary. Each noise signal should have half of the calculated variance as, when combined, this will result in a signal with our total calculated variance. Another way to think about this is in terms of power. When we have zero mean, power and variance are equal. Therefore, if we have a complex signal of a certain power, this will comprise two independent signals each with half of the stated power.\n",
    "\n",
    "We can generate the real and imaginary components of the noise signal using *np.random.normal()*. This function takes three arguments: mean, standard deviation and size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_noise = math.sqrt(var_noise/2)     # Standard deviation of 1-D noise signal\n",
    "mean_noise = 0\n",
    "len_noise = len(signal)\n",
    "\n",
    "noise_I = np.random.normal(mean_noise,std_noise,len_noise)\n",
    "noise_Q = np.random.normal(mean_noise,std_noise,len_noise)\n",
    "\n",
    "noise = noise_I + 1j*noise_Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the cell below, we can verify that our noise signal has a variance close to that which we calculated given a specific SNR. Again, the more bits that this is calculated over, the more accurate this result will be. If one of the shorter codes is being employed, the measured variance can differ from the desired variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Expected Noise Variance:',var_noise)\n",
    "print('Measured Noise Variance:',np.var(noise))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below plots the noise signal. Try varying the desired SNR above and re-running the cells. We can see how changing the SNR value results in a different amplitude range. The larger the SNR the smaller the noise amplitudes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hf.plot_samples('Noise Signal',\n",
    "                [range(len(noise_I)),range(len(noise_Q))],\n",
    "                [noise_I,noise_Q])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also inspect the PDF of our noise. The blue histogram shows our measured PDF whereas the red line shows the ideal Gaussian PDF given our desired standard deviation. The parameter *b* can be varied to change the number of bins, or the resolution along the x-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hf.plot_hist('',np.concatenate((noise_I,noise_Q)),std_noise,b=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. Add Noise <a class=\"anchor\" id=\"add_noise\"></a>\n",
    "To add the generated noise to our signal, we simply add the two arrays together. The cell below performs this addition and plots the encoded noisy signal's samples  as well as a constellation plot. Inspecting these plots, it can be seen as the SNR decreases the ability to determine which symbol a sample belongs to becomes more difficult. This is how errors are introduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal_with_noise = signal + noise\n",
    "\n",
    "hf.plot_samples('Encoded Noisy Signal',\n",
    "                [range(len(signal_with_noise)),range(len(signal_with_noise))],\n",
    "                [np.real(signal_with_noise),np.imag(signal_with_noise)])\n",
    "\n",
    "hf.plot_constellation(signal_with_noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can verify that the correct level of noise has been added to the signal by measuring the SNR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "power_signal = np.var(signal)\n",
    "power_noise = np.mean(np.abs(signal - signal_with_noise)**2)\n",
    "snr_db = 10*np.log10(power_signal/ power_noise)\n",
    "print('Desired SNR: ', SNR)\n",
    "print('Measured SNR: ', snr_db)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 4. Log Likelihood Ratio (Soft Demodulation) <a class=\"anchor\" id=\"llr\"></a>\n",
    "Once a signal has been received, it must be baseband demodulated. This is the opposite process to baseband modulation. We receive a signal with noise and must determine which symbol a sample belongs to so that we can convert that symbol back into the original serial sequence. If there is no or little noise in a channel, then this process can be relatively simple. Boundaries around the symbols can be drawn up and if a sample falls within a boundary's region, then it must belong to that symbol. An example of this is shown in Figure 3 below.\n",
    "\n",
    "<a class=\"anchor\" id=\"fig-3\"></a>\n",
    "<center><figure>\n",
    "<img src='./images/scaled_symbol_mapping_boundaries.svg' width='800'/>\n",
    "    <figcaption><b>Figure 3: Gray coded symbol mapping constellation for 16-QAM with decision boundaries.</b></figcaption>\n",
    "</figure></center>\n",
    "\n",
    "This straight conversion from a sample back to a symbol can be referred to as a **hard** decision and results in hard bits. There is no nuance in the decision made. A sample could be right on the boundary between symbols and, after baseband demodulation, the resulting bits would have the same perceived certainty that the bits resulting from a sample that aligns perfectly with a symbol has. Therefore, making hard decisions can lead to errors.\n",
    "\n",
    "When decoding data and performing error correction, it is beneficial to instead make **soft** decisions during this process. Soft decisions return soft bits which convey the confidence that a bit is either a 1 or a 0. In doing this, the decoder is able to form a better estimate of the original data which can result in fewer errors than if hard decisions were employed. \n",
    "\n",
    "The name Soft Decision Forward Error Correction (SD-FEC) of the integrated blocks found in some Xilinx RFSoC devices means it should be of no surprise that soft bits are the form accepted by the decoder. Soft bits themselves can come in different formats. The particular form of soft bit accepted by the SD-FEC core is **Log-Likelihood Ratio (LLR)**. The LLR is calculated for each bit in a symbol and is defined as: \n",
    "\n",
    "$$\n",
    "    \\mathit{LLR(b)} = ln\\left(\\frac{Pr(b=1)}{Pr(b=0)}\\right)\\\\\n",
    "$$\n",
    "This means that positive values (including 0) are interpreted as a hard binary 1 and negative values are interpreted as a hard binary 0. The magnitude of these values indicates the certainty, with a larger magnitude equating to a higher certainty and vice versa.\n",
    "\n",
    "This calculation can be reduced to a series of if-statements for easy implementation. Krishna Sankar has written an excellent [blog](http://www.dsplog.com/2009/07/05/softbit-16qam/) on how soft bits can be calculated when demapping 16-QAM signals. For our 16-QAM mapping, the resulting conditional statements are:\n",
    "\n",
    "$$\n",
    "\\begin{gather}\n",
    "  llr(b3) = \n",
    "  \\begin{cases}\n",
    "      2(y_{re}+1) &\\;\\;\\;: \\text{if}\\; y_{re} < -2\\\\\n",
    "      \\;\\;\\;\\;y_{re} &\\;\\;\\; : \\text{if}\\; -2 \\leq y_{re} < 2\\\\\n",
    "      2(y_{re}-1) &\\;\\;\\; : \\text{if}\\; y_{re} > 2\n",
    "  \\end{cases}\n",
    "\\end{gather}\n",
    "$$\n",
    "<br>\n",
    "<br>\n",
    "$$\n",
    "\\begin{gather}\n",
    "  llr(b2) = \\;\\; -|y_{re}|+2 &&\\; \\forall \\; y_{re}  \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\n",
    "\\end{gather}\n",
    "$$\n",
    "<br>\n",
    "<br>\n",
    "$$\n",
    "\\begin{gather}\n",
    "  llr(b1) = \n",
    "  \\begin{cases}\n",
    "      -2(y_{im}+1) & : \\text{if}\\; y_{im} < -2\\\\\n",
    "      \\;\\;\\;\\;-y_{im} & : \\text{if}\\; -2 \\leq y_{im} < 2\\\\\n",
    "      -2(y_{im}-1) & : \\text{if}\\; y_{im} > 2\n",
    "  \\end{cases}\n",
    "\\end{gather}\n",
    "$$\n",
    "<br>\n",
    "<br>\n",
    "$$\n",
    "\\begin{gather}\n",
    "  llr(b0) = \\;\\; -|y_{im}|+2 &&\\; \\forall  \\; y_{im} \\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\\;\n",
    "\\end{gather}\n",
    "$$\n",
    "\n",
    "This is easily implemented in Python. The cell below iterates through each sample in the noisy signal and returns 4 LLR values, one for each bit in the symbol, for each sample in the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "llrs = []\n",
    "for y in signal_with_noise:\n",
    "    # Bit 3\n",
    "    if y.real < -2:\n",
    "        b3 = 2*(y.real+1)\n",
    "    elif y.real >= -2 and y.real < 2:\n",
    "        b3 = y.real\n",
    "    elif y.real > 2:\n",
    "        b3 = 2*(y.real-1)\n",
    "        \n",
    "    # Bit 2\n",
    "    b2 = -abs(y.real)+2\n",
    "    \n",
    "    # Bit 1\n",
    "    if y.imag < -2:\n",
    "        b1 = -2*(y.imag+1)\n",
    "    elif y.imag >= -2 and y.imag < 2:\n",
    "        b1 = -y.imag\n",
    "    elif y.imag > 2:\n",
    "        b1 = -2*(y.imag-1) \n",
    "        \n",
    "    # Bit 0\n",
    "    b0 = -abs(y.imag)+2\n",
    "    \n",
    "    llrs.append(b3)\n",
    "    llrs.append(b2)\n",
    "    llrs.append(b1)\n",
    "    llrs.append(b0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visually examine how the soft decisions have performed by comparing the first 80 LLR values with the first 80 bits in the uncorrupted serial data. Executing the cell below will plot these two signals. The blue plot shows the LLR values and the red the serial data before noise was added. Remember that positive and zero LLR values indicate a 1 and negative LLR values a 0. The LLR magnitude indicates the certainty so if any errors are observed here, the LLR value should have a small magnitude. Try changing the SNR and running the previous cells again to see how this plot is affected by noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hf.plot_samples('Comparing LLR Values to Actual Bits',\n",
    "                [range(80),range(80)],\n",
    "                [llrs[0:80], rx_bits[0:80]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, store the LLR values so that they can be used in the following notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%store llrs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Conclusion <a class=\"anchor\" id=\"nb3_conclusion\"></a>\n",
    "This notebook has shown a simplified radio communications pipeline, where encoded serial data has been mapped into symbols, subjected to a noise channel, and demapped using soft decisions into soft bits. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "[⬅️ Previous Notebook](02_fec_encoding.ipynb) || [Next Notebook 🚀](04_fec_decoding.ipynb)\n",
    "\n",
    "Copyright © 2023 Strathclyde Academic Media\n",
    "\n",
    "---\n",
    "---"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
